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ABSTRACT 

We perform a systematic study of how the inhomogeneities in the intergalactic medium 
(IGM) affect the observability of Lya emitters (LAEs) around the epoch of reioniza- 
tion. We focus on the IGM close to the galaxies as the detailed ionization distribution 
and velocity fields of this region could significantly influence the scattering of Lya 
photons off neutral H atoms as they traverse the IGM after escaping from the galaxy. 
f£h We simulate the surface brightness (SB) maps and spectra of more than 100 LAEs 

Qh at z = 7.7 as seen by an observer at z = 0. To achieve this, we extract the source 

q properties of galaxies and their surrounding IGM from cosmological simulations of 

box sizes 5 — 30 h~ 1 Mpc and then follow the coupled radiative transfer of ionizing 
^ and Lya radiation through the IGM using CRASHa. We find that the simulated 

SB profiles are extended and their detailed structure is affected by inhomogeneities 
in the IGM, especially at high neutral fractions. The detectability of LAEs and the 
fraction of the flux observed depend heavily on the shape of the SB profile and the 
SB threshold (SB^) of the observational campaign. Only ultradeep observations (e.g. 
^j- SB^ ~ 10 -23 erg s _1 cm -2 arcsec -2 ) would be able to obtain the true underlying 

mass-luminosity relation and luminosity functions of LAEs. The details of our results 
\/~) depend on whether Lya photons are significantly shifted in the galaxy to longer wave- 

C\( lengths, the mean ionization fraction in the IGM and the clustering of ionizing sources. 

These effects can lead to an easier escape of Lya photons with less scattering in the 
IGM and a concentrated SB profile, similar to the one of a point source. Finally, we 
show that the SB profiles are steeper at high ionization fraction for the same LAE 
^— H sample which can potentially be observed from the stacked profile of a large number 

of LAEs. 
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1 INTRODUCTION 



Partridge &; Peebles| ( |1967| predicted that high redshift 
galaxies would emit copious amounts of Lya photons (^10 
per cent of the galaxy luminosity) due to the reprocessing 
of ionizing photons from the young massive stars by the in- 
terstellar medium (ISM). This stimulated a large number of 
studies aimed at observing the high-redshift galaxies which 
mostly lead to non-detections (e.g.|Davis &; Wilkinson|1974) 
|Partridge|1974| |Djorgovski & Thompson||1992| and the ref- 
erences there in). Models were made to reason the null de- 
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tections such as absorption of Lya radiation by dust (e.g. 
| Spitzer|1978l|Meier fe Terlevich|198T]|Hartmann et al.|1988| 
|Charlot &; Fall|1993| and lower number of massive stars due 
to ageing of stellar populations ( Valls-Gabaud 1993). 

The first detections of high redshift Lya emitting galax- 
ies (which we also refer to as 'Lya Emitters', or LAEs) were 
reported in the mid-1990s (e^. |Hu fe McMahon|1996| ). Since 
then a large number of LAEqjhave been observed using both 
narrow band photometry and spectroscopy up to z < 7 (e.g. 



1 The term 'LAE' often refers to galaxies that have been selected 
on having a 'strong' Lya emission line. In this paper we use the 
term LAE for all Lyo; emitting galaxies. Our classification thus 
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Rhoads fe Malhotra||200Tl |Dawson et al.||2004| |Rhoads et l 
al.|2004||Taniguchi et al.|2005||Iye et al.|2006||Nilsson et aT] 
2007[|2011||Ouchi et aI.|20091|2010||Wang et al.|2009||Guaita 



et al.|2Q10[|Lehnert et al.|2010l |Wolf et al.|2010| ICassata et' 
al.||2011||Gronwall et al.||20111 jkashikawa et al.|2011| |Pen- 
tericci et al.||2011| |Ono et al.||2012| |Schenker et al.||2012^ 7 

In addition to being one of the main methods of detecting 
high-redshift galaxies to study the evolution of galaxy prop- 
erties with redshift, LAEs can also be used to probe the 
epoch of reionization (EoR) 



(e.g. 



Miralda-Escude & Rees 



1998; |Malhotra fe Rho ads 2004; |McQuinn et al.||2007| ) and 



to constrain dark energy properties using baryonic acoustic 



oscillations in the LAE clustering signal (e.g. Eisenstein et 
al.| |2005| |Hill et al.||2008| ). Since this paper is motivated by 
the use of LAEs for the detection of EoR, we will discuss 
this in more detail. 

Due to its resonant nature, the Lya line is scattered by 
even very small amounts of neutral hydrogen along its path. 
The optical depth r due to scattering away from the line-of- 
sight to the observer (at z = 0) for a LAE at redshift z s is 
given by 



6.02 x lu x m ( ' 



3/2 



(i) 



for photons entering the intergalactic medium (IGM) with 
wavelength A < A a = 1215.67 A (Lya rest wavelength; e.g. 
Gunn & Peterson 1965; Barkan a" &; Loeb||2001| and by 



r 2.9 xm 



Av 



600 km s 



1 + Zs 

10 



3/2 



(2) 



A a (l + Av/c) for 



for photons entering the IGM with A 
an equivalent velocity shift of Av (e.g. Miralda-Escude & 
|Rees||1998| |Dijkstra fc WyithepOlO) . Here, x m is the neu- 
tral fraction of the IGM, c is the speed of light in vacuum 
and the cosmology used to derive the above equations and 
for the rest of the work in this paper is the standard model 
- Q ,m=0.3, Q ,a=0.7, Q , 6=0.04, h=0.7, n=l and cr 8 =0.9. 
The effect of scattering on the shape of the intrinsic spec- 
trum can be used to estimate the ionization fraction of the 
IGM around the source and thus study the EoR. 

LAEs luminosity functions (LFs), number density and 
clustering are the main methods used to study the EoR his- 
tory. It is expected for example that an increment in the 
neutral fraction of the universe leads to a reduction of the 
observed number of sources, thus supressing the LFs towards 
higher redshifts 



(e.g. 



Haiman & Spaans 1999| . Following 



this idea, |Ouchi et al. (2010) used the drop in LAEs LF 
at z = 5.7 — 6.6 to constrain xm < 0.2 ± 0.2 at z = 6.6. 
A less stringent limit of xm ^$0.5 has been found at the 
same redshift by |Ouchi etaT] ( |2Q10| ) using measurements of 
LAE clustering. These are marginally consistent with the 
estimate by Taniguchi et al.| ( |2005| that at least 20-50 per 
cent of the volume of the universe needs to be ionized to 
account for the observed number of LAEs at z = 6.5. The 
use of number density (e.g. Malhotra & Rhoads 2006) and 
clustering p roperties (e.g. |Furlanetto, Zal darriaga, & Hern- 
quist||2006l |McQuinn et al.||2007|) of LAEs have also been 



includes ultraviolet (UV) selected galaxies for which follow-up 
spectroscopy revealed the presence of a strong Lyo; emission line. 



extensively discussed in the literature as tools to probe the 
EoR. 

A large number of authors have modelled LAEs at dif- 
ferent redshifts using analytic (e.g. [Ha iman 2002 ; Kobayashi, 



Kamaya, &; Yonehara 2006; Dijkstra, Wyithe, &; Haiman 



2007| [K obayashi, Totani, & Nagashima 2007; Sa mui, Sri- 



anand, & Subramanian 2009 Tilvi et al. 2009) and semi- 



numeric (e.g. |Le Delliou et ai.||2005| |2006| |Furlanetto, ZaT- 



darriaga, & Hernquist 2006; Tasitsiomi 2006; McQuinn et 



aLl|2007l pev et al. 2008; Dayal et al.||2009| [Zheng et al. 



2010; Dijkstra, Mesin ger, fe Wyithe| [2011] |Yajima et al. 



2011 ) methods. It has been shown that properties like dust 
content (e.g. Kobayashi, Totani, & Nagashima 2007; Dayal 



et al. 2009) and gas inflows (e.g. |Santos 2004; Dijkstra, 



Lidz, & Wyithe 2007) reduce the observability of LAEs, 
whereas top-heavy initial mass function (IMF; e.g. |Malho- 



tra fe Rhoads||2002| |Dijkstra fe Wyithe||2007|), clustering of 



ionizing sources, 'patchy' reionization (e.g. |Wyithe &; Loeb 



2005; Furlanetto, Zaldarriaga, & Hernquist 2006; McQuinn 
et al.|2007||Iliev et al.|2008||Mesinger fe Furlanetto|2008a|bf 
Dijkstra, Mesinge r, &; Wyithe||2011[) and gas outflows (e.g. 
Santos|2004||Dijkstra fe Wyithe|2010| ) improve their observ- 



ability at each redshift. Most works simulate the scattering 
of Lya photons away from the line of sight by suppressing 
the intrinsic spectrum of an LAE by the optical depth calcu- 
lated along the line of sight (see equations [I] and [5]) . Only a 
few studies perform a full 3D RT (RT) of Lya photons (e.g. 
|Cantalupo et al. 2005; Tasitsiomi 2006; Laurs en~&; Sommer^] 
|Larsen|[2007| |Laursen, Sommer-Larsen, &; Andersen 2009 



Kollmeier et al.|2010| |Zheng et al.|2010| |Barnes et al.|20TT 
|Yaj 



ima et 



al.|2011 ) 



In this work, we look at the structure in the IGM at 
z = 7.7, where the current constraints point towards non- 
negligible neutral fraction in the volume averaged IGM, 
specifically focusing on the effects of the gas close to the ob- 
ject. This region is important because the inhomogeneities 
in the gas density play a crucial role in determining the 
size of the ionized region, especially at very low ionization 
fractions. The residual neutral hydrogen close to the source, 
together with the velocity field in this same region strongly 
affects the scattering pattern of the Lya photons and their 
escape towards an observer. Here we investigate these issues 
for the first time in detail for a large (i.e. > 100) sample of 
objects. LAEs have been detected up to z ~ 7 (e.g. |Iye et| 
al.|2006 Ono et al.|2012 ). Because observations in the near 
future are aiming for z - 7.7 (e.g. |Tilvi et al.|2010|pement[ 



et al.|20l2| , in this study we have chosen to concentrate on 
objects at z = 7.7. 



Zheng et al. (2010) previously simulated a large sample 



of LAEs with Lya RT at moderate spatial resolution in a 
highly ionized IGM at z — 5.7. Similarly, Laursen, Sommer- 



Larsen, & Razoumov (2011) studied scattering in the post- 



reionization IGM using hydrodynamical simulations that 
had much higher spatial resolution for RT. Their work is 
complementary to our study. 

In this work we simulate LAEs using galaxies from 
hydrodynamical cosmological simulations and RT using 
CRASHa ( jPierleoni, Maselli, fe Ciardi|2009| ). We make sur- 
face brightness (SB) maps of the objects and study the effect 
of IGM structure on the SB profiles. In the following Sec- 
tion, we describe the cosmological simulations and Section [3] 
briefly introduces the RT code CRASHa. Section [4] explains 



© 2012 RAS, MNRAS 000, [Tp2] 



Effect of IGM on LAE detections during EoR 3 



Model 


L (/i~ 1 Mpc) 


number of particles 


^DM (M ) 


^gas 


(M ) 


r\ (h x kpc) 


L05 


5 


2 x 320 3 


3.93 x 10 5 


6.04 


x 10 4 


0.78 


L10 


10 


2 x 320 3 


3.14 x 10 6 


4.83 


x 10 5 


1.56 


L20 


20 


2 x 320 3 


2.52 x 10 7 


3.87 


x 10 6 


3.13 


L30 


30 


2 x 320 3 


8.49 x 10 7 


1.30 


x 10 7 


4.69 



Table 1. Simulation properties. From let to right: model name; comoving box size, L; total number of particles (DM and gas); mass of 
DM particles, ttzdm; mass of gas particles, m gas ; softening length, 77. 



the pipeline used to simulate a statistically significant sam- 
ple of LAEs. The results are described in Section [5] and the 
parameter study to understand the uncertainties in our work 
is shown in Section [6] In Section [7] we discuss the effect of 
IGM structure on the methods used in estimating the mean 
ionization fraction of the universe. Finally the conclusion are 
drawn in Section [8] 



2 SIMULATIONS OF GALAXY FORMATION 

To study high redshift galaxies and their surrounding IGM, 
we need simulations with a large box size to provide a sta- 
tistical sample of halos in a wide mass range in different 
environments. At the same time, high resolution to resolve 
the halo and the structure in the surrounding IGM is re- 
quired. A compromise was achieved by using medium range 
box sizes of 5-30 /i _1 Mpc comovine;. The simulations used 



THE RADIATIVE TRANSFER CODE 
CRASHa 



in this paper are described in Maio et al. (2010), although 



additional ones have been run for different box sizes. 

The simul ations were run using the TREE-PM SPH 
code Gadget-2 |Springel 2005) modified to include primor- 
dial Hydrogen, Helium and Deuterium based chemistry [e~, 
H, H + , He, He + , He ++ , H 2 , H+, H~, HeH + , D, D + , HD] 



( Yoshida et al.|2003||Maio et al.|2007[), s tellar evolution and 
metal pollution (Tornatore et al. 2007| ) and fine structure 



metal transition cooling (O, C + , Si + , Fe + ) at T < 10 4 K 



(Maio et al. 2007 2009). Metals produced by AGB stars 



and supernovae (SNII, SNIa) are spread by supernovae/wind 
feedback. The IMF is chosen according to the metallicity of 
the stellar particles, Z. We assume that a transition from 
metal- free/ very metal-poor Population III stars to metal- 
enriched, more standard Population II/I stars takes place 
whenever the gas reaches a critical metallicity Z cr it, which 
is determined by the cooling and fragmentation properties 
of the gas. Z cr it quoted by different authors range between 
10~ 6 to 10~ 3 - 5 Z Q (e.g. 



2003 



Bromm & 



Schneider et al. 
Loeb|2003| [Schneider et al.|2006| |Santoro"fe Shuil|2006| ). In 
these simulations we choose Z cr it = 10~ Z©. Below Z cr it 
we assume a Salpeter IMF in the mass range [100,500] M©, 
while for metallicities above Z cr i t a Salpeter IMF in the mass 
range [0.1,100] M© is used. 

The simulations analyzed in this paper have box sizes 
of L = 5, 10, 20 and 30 /i _1 Mpc comoving with 320 3 parti- 
cles each in dark matter and gas. The mass of dark matter 
particles is 3.93 x 10 5 M© (L/5 /i _1 Mpc) 3 and the gas par- 
ticle mass is 6.04 x 10 4 M© (L/5 /i^Mpc) 3 . The comoving 
softening length is 0.78 /i -1 kpc (L/5 h~ 1 Mpc) 3 which is 
~ 1/20 of the mean inter-particle distance. The details of 
the simulations are given in Table [T] 



CRASHa flPierleoni, Maselli, fe Ciardi|2009| ) is the first RT 
code for cosmological application where the propagation of 
Lya and ionizing photons are coupled (refer to Yajima et] 



al. 2011 for a similar code). The ionizing part is based on 
CRASH JcTardi et al.|2001|jMaselli, Ferrara, fe Ciardi|2003 



Maselli, Ciardi, & Kanekar 2009) which is a 3D ray-tracing 



grid-based RT code using Monte-Carlo (MC) techniques 
to sample the probability distribution of several parame- 
ters like spectrum of sources, emission direction and optical 
depth. The ionizing radiation is propagated through an ar- 
bitrary static H/He gas density field. Both radiation from 
point sources and diffuse radiation from ultraviolet back- 
ground or recombination of H/He gas can contribute to the 
ionizing flux. 

The total ionizing photon count, E s [phot], from each 
sources with photon rate iVion,s [phot s _1 ] over the whole 
time tsim of the RT simulation, is distributed among N p pho- 
ton packets containing ionizing photons of different frequen- 
cies depending on the source/background spectrum. The 
emission of photon packets from each source happens at 
equally spaced time intervals dt — t s i m /N p . The direction 
of emission is assigned by the MC sampling of the angular 
characteristic of the source. As the packets are propagated 
through each cell z, photons are absorbed according to the 
cell optical depth rl calculated combining the contributions 
from HI, Hel and Hell. The probability that a photon in a 
packet is absorbed is calculated using P{r l c ) — 1 - e T °. This 
changes the temperature and ionization conditions of the 
cell and the photon distribution in a packet. The trajectory 
of a packet is followed till all the photons are absorbed or in 
case of non-periodic boundary conditions, till it reaches the 
edge of the box. 

CRASH is modified to follow the time evolution of 
Lya photons through an evolving ionization configuration 
of gas. A statistical approach to the Lya RT is adopted us- 
ing pre-compiled tables from MCLya ( |Verhamme, Schaerer, 
& Maselli 2006). Lya radiation is emitted by both point 



sources and recombination of ionized gas. Just like in the ion- 
izing case, the Lya photons produced by each point source s 
with photon rate Ni^ ya ^ over the entire simulation time t s i m 
is divided between Lya photon packets. N p ^ ya Lya photon 
packets are emitted by each source s at regular time inter- 
vals dti,y a = tsim/^V e m,i in random directions. Af em ,i is the 
number of times the Lya photons are propagated during the 



2 Each galaxy within the hydro-dynamical simulations which con- 
tains at least a star forming particle is represented within the RT 
simulation as a point source. 
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Model 


M B m (x 


10 10 M ) 


M* (x 10 7 M ) 


M gas (x 10 8 M ) 


Z (x 10 


- 2 z ) 


L05 


0.12 


- 2.58 


6.61 - 256 


0.03 - 3.63 


1.20 - 


6.87 


L10 


1.33 


- 9.70 


9.06 - 108 


15 - 110 


2.02 - 


4.80 


L20 


1.84 


- 8.57 


4.59 - 64 


26 - 110 


0.76 - 


3.52 


L30 


3.58 


- 16.6 


3.54 - 134 


47.2 - 213 


0.34 - 


2.85 



Table 2. Properties of haloes at 2=7.7. From left to right: name of the model; range of DM halo mass, Mdm; stellar mass, M*; gas 
mass, Mga S ; gas metallicity, Z. 



RT simulation. When a packet enters a cell, depending on 
the conditions in the cell, the photons are either absorbed 
or scattered. In case of scattering, a new wavelength and 
the time of escape from this cell is obtained from the pre- 
calculated tables. The typical time of escape from a cell of 
size 1 kpc, temperature 10 K and number density 0.01 cm -3 
giving an optical depth of ~ 6 x 10 7 is about 3 x 10 5 years, 
(e.g. |Adams|1975l |Dijkstra fc Loeb|2008t Photons from the 
recombination in the gas are added to the packets. 

Dust is an important factor which determines Lya ra- 
diation transport (e.g. Neufeld 1991). Dust optical depth is 
calculated as a fraction f T of the Lya optical depth. f T is de- 
termined as f T = m p / dust x in/dust x cr du st x p ce ii x d ce ii where 
proton-to-dust mass ratio m p / dust — 5.0 x 10~ 8 , gas-to-dust 
ratio /n/dust — 5 x 10 ~ 3 (Verhamme, Schaerer, & Maselli 



2006 ) , dust absorption cross section adu 



2 x 7T x r dust for 
a dust grain of radius r dust = 2.0 x 10 -(3 cm. p ce u is the gas 
density in the cell and d ce \\ is the distance traveled within 
the cell by the photon. In high density regions, where Lya 
photons have a high probability of scattering multiple times 
within the cell volume, the probability of dust absorption 
gets high, with the consequence that the effects due to dust 
scattering on Lya RT are negligible. For example, assuming 
that the probability of scattering of Lya photons by dust is 
equal to that of absorption, in 20 photon-dust interactions, 
only 1 in 10 6 photons would continue to scatter without be- 
ing absorbed. Thus in high density regions, we neglect dust 
scattering. 

For more details, we refer the reader to the original 
papers. 



4 SIMULATING LYa EMITTERS 

Our plan is to study how the observability of LAEs is af- 
fected by transmission through the IGM. For this, we calcu- 
late Lya SB profiles of a large number of simulated LAEs at 
z — 7.7 covering a wide range of dark matter halo masses, 
thus sampling an equally wide range in IGM environments. 
The redshift was chosen as observational efforts are under- 
way to investigate LAEs at z 



7.7 (e.g. |Tilvi et al.||2010 



Clement et al.|2012| ). As at this relatively high z the IGM is 
expected to be still substantially neutral, a full 3D RT ap- 
proach is necessary to investigate the observability of LAEs. 

The method adopted to simulate LAEs is to extract a 
cube around each dark matter halo from the snapshot of the 
simulation of galaxy formation, grid the density and veloc- 
ity fields, get the spectrum produced by its stellar popula- 



ionization fields and run the RT simulation to obtain Lya 
photons escaping the simulation cube. We can use the de- 
tails of the photons exiting each of the six faces of the cube 
to make SB maps. Each step is explained in detail in the 
following Sections. 

4.1 Extracting Halos and Gridding 

A snapshot of the hydrodynamical simulation is taken and 
the properties of the dark matter haloes are obtained us- 
ing the Friends-of- Friends (FoF) routine ( [Davis et al.|1 985). 
The centre-of-mass (CoM) of the dark matter (DM) halo, 
the DM, gas and stellar mass, and mean gas metallicity are 
used in the rest of the analysis. The mean gas metallicity 
is defined as the fraction of gas mass in metals. A sum- 
mary of the DM halo properties for the four different simu- 
lations used in this paper are given in Table [2] In our study, 
we choose a subset of the resolved haloes (i.e. with mass 
M DM > 100 771dm 5 e.g. |Trenti et aC||2010| ) from L05-30 to 



get a fair sampling of the DM mass range 1O 9_11 M0. This 
results in ~ 130 haloes. 

In our simulations, we resolve the structure in the IGM, 
but we do not properly resolve the ISM of the individual 
galaxies. To mimic the absorption/scattering of radiation 
within the ISM, we use the commonly adopted prescription 
based on the escape fraction of ionizing photons, /esc, ion, ism, 
and the escape fraction of Lya photons, / eS c,Lya,iSM (see 
Sec. 6.1 ). Thus we need to exclude from our RT calculations 



tion from STARBURST99 ( [Leitherer et al.|1999) using the 
halo properties, input all the details to CRASHa along with 
the default values (refer to Section 4.3) for temperature and 



those cells which represent the ISM, to avoid accounting for 
their effect twice. We obtain this by removing all cells with a 
density larger than 0.03 cm -3 . Even values up to 0.05 cm -3 
are acceptable, but might show slight dependence on the 
gridded density resolution especially for high resolution grid- 
ding in the case of small haloes (see discussion in the next 
paragraph). For the case of 0.03 cm -3 , the effect of density 
resolution is negligible in the range of simulations we deal 
with in this work. Values < 0.03 cm -3 give a similar ion- 
ization structure but lead to loss of some of the IGM gas. 
As a precaution to avoid removal of gas from high density 
regions not associated with the specific dark matter halo, 
we restrict the removal of cells to a distance of ~ 0.7 x r2oo 
from the source. The precise choice of this radial distance 
does not affect our main results, as long as it is less than the 
nearest massive dark matter halo in the cube. 

A cube around the CoM of the halo is extracted with a 
side of 35 x r2oo, where r2oo is the radius in comoving units 
at which the mean DM density inside the sphere is 200 times 
the critical density of the universe at that redshift. This is 
the smallest box containing the HII region produced by the 
galaxy at its equilibrium, so that we get the highest possible 
spatial resolution for the gridding of density and velocity 
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fields are gridded as input for CRASHa. For gridding, we use 
GadgettoGrid ( Pakmor|2010 ). The code distributes the par- 
ticle mass and peculiar velocities using an SPH kernel to 64 
neighboring particles and then grids the fields. The default 
grid size used is 256 3 , which is set by memory and runtime 
constraints of CRASHa, although different grid sizes have 
been used for testing purposes. An important caveat is that, 
because the grid resolution is a factor of the virial radius of 
the halo, this results in a lower physical resolution for higher 
mass objects. 

Figure [I] shows a slice through the simulation box. The 
top panel refers to the gas number density in the box with 
the to-be-extracted cube around a halo marked by a solid 
red line and the CoM of the halo marked by a green dot. The 
extracted density field is gridded and shown in the bottom 
panel. Note that the gas density in the IGM is not uniform 
and the gridded density field has a range in values spanning 
many orders of magnitude. 



Gas Number Density fcm 3 ] 
-5 -4 -3 -2 -1 




500 1000 1500 2000 
X [kpc/h] 

Figure 1. Example of the extraction of a cube around a DM 
halo and the gridding of the gas density and velocity fields for 
CRASHai. Top: The gas distribution within a thin slice of the 
simulation box L05 centred on the CoM (shown as a green dot) 
of the most massive halo at z = 7.7 in comoving units of length. 
The comoving thickness of the slice is 10 h~ 1 kpc. The cube to be 
extracted around the object for the RT simulation is shown with 
a solid red square. Bottom: The gridded gas number density field 
(in logarithmic scale) corresponding to the cube in the top panel. 



fields for our RT simulations. The exact value of the cube 
dimension though does not affect our main results. A cube 
size in the range 25 — 50 x r2oo gives similar results, but the 
lower bound is too close to the edges of the ionized region, 
while larger boxes would have lower spatial resolution at 
fixed grid size, which is undesirable. 

Once the cube is extracted, the gas density and velocity 



4.2 Luminosity of Stellar Sources 

As mentioned in the previous Section, we obtain the halo 
properties using the FoF and use the associated IMF, stel- 
lar mass and metallicity to calculate the corresponding ion- 
izing and Lya photon rate. To this purpose, we use STAR- 
BURST99 (Leitherer et al. 1999), which is a population syn- 
thesis code that provides the spectrum of a stellar popula- 
tion given an IMF, stellar mass, metallicity and age. 

The lowest metallicity available in STARBURST99 is 
Z = 0.0004 - 0.02 Z for Padova tracks. We choose the 
Padova original stellar evolution tracks instead of Padova 
AGB tracks as the Hubble age at z = 7.7 is less than 1 
Gyr, allowing us to ignore the contribution of AGB stars to 
the spectrum. The average ionizing flux for the same stel- 
lar mass is mildly higher for lower m etallicities, i.e. ~ 17% 
higher for Z = 3 x 10~ 3 Z (refer to |Schaerer||2003| for low 
metallicity values) . Since our metallicities are above Z CT [ t , we 
use a Salpeter IMF in the mass range [0.1,100] M in the 
STARBURST99 models, consistently with our cosmological 
simulations. The spectrum is scalable with respect to the 
stellar mass of the halo. Therefore, we choose M* = 10 6 M 
as the default stellar mass for the STARBURST99 model 
and normalize the photon count to this stellar mass. We 
also choose instantaneous star formation mode for comput- 
ing the spectrum. The mean ionizing photon rate Ni on was 
obtained by averaging the number of photons emitted over 
t = 2 X 10 years in the wavelength range [91, 912] A. The 
time t was chosen as the cumulative number of ionizing pho- 
tons reaches convergency. Thus, we obtain a normalized ion- 
izing photon rate iY ion = 7.77x 10 50 (M*/10 6 M ) phot s" 1 . 
Only a fraction of these ionizing photons reaches the IGM, 
the rest is converted to Lya photons in the ISM. Therefore, 
the ionizing photon rate reaching the IGM, iVJfn, is: 



N e 

1 v ir 



/esc, ion, ISM X N[ c 



(3) 



where /esc, ion, ism is the escape fraction of ionizing photons 
from the ISM. Our reference value is /esc, ion, ism = 0.02 



(Gnedin, Kravtsov, & Chen 2008). In Section 6.1 we dis- 



cuss more extensively this choice. 

The calculation of the Lya photon rate iVLyc* is more 
complicated. There are three different components consid- 
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Figure 2. Sketch of the method used to calculate the SB. The star represents the pixel with the stellar source. The dotted lines represent 
the radial direction connecting each pixel on the edge of the box to the source. The solid lines with arrows represent the directions of 
photon packets within such pixels. Dashed lines connect the pixels to the observer at z = 0. 



ered in our work - stellar continuum Nhya,*, nebular contin- 
uum A^Lya, neb and recombination in the ISM ALya, ism: 



A;Lya — N-Lya,* + A^Lyc*, neb + ^Lya, ISM- 



(4) 



Nhya,* and A^ya, neb are calculated from the spectrum 
by averaging over t at 1215.67A with a line width of 2 A, 



which is the typical intrinsic line width of LAEs ( Partridge 
fe Peebles! [1967[ . In STARBURST99 spectra, the nebular 



continuum is calculated by converting all the ionizing pho- 
tons into nebular lines except for Lya. Therefore, Nh ya , neb 
gives the upper limit of the contribution from other nebular 
lines to this wavelength range. 

A^Lya, ism is instead computed from the ionizing photon 
rate as (also see e.g. |Schaerer|2003l ): 



N] 



Lya, ISM 



0.68 x N ion x (1 - /« 



esc, ion 



ism). 



esc, ion, ISM - 
,48 



Assuming our reference value of /< 
of Lya photons are N-L ya ,* = 6.94 x 10 
phot s~\ A Lya , n eb = 2.83 x 10 47 (M*/10 6 
and N~Lya, ism = 5.07 x 10 50 (M*/10 6 M ) 



(5) 

= 0.02, the rate 
(M*/10 6 M ) 
M ) phot s" 1 
phot s _1 (for 

— K 1 

ya. 



an ionized gas temperature of 10 4 K). Thus, Nt, ya = 5.15 x 
10 50 (M*/10 6 M ) phot s -1 and our total Lya photon rate 
is completely dominated by photons that were emitted as 
recombination radiation in the ISM. 

Finally, the Nf^ a which escapes the galaxy into the 
IGM after absorption by dust can be denned as: 



A^Lya — /esc, Lya, ISM X A^Lyo 



(6) 



where /esc, Lya, ism is the escape fraction of Lya after absorp- 
tion by dust in the galaxy. Our fiducial value is /esc, Lya, ism = 
0.3 flDayal et aL"1|2009| ). In Section |rH we discuss more ex- 
tensively this choice. N^ y c a can also be converted into a lu- 
minosity, ^Lya, by convolving with the spectrum of Lya 



photons. 

Our value of N- 1C 
bers quoted by e.g. iDayal et al. (2009) and Zheng et al. 



is a factor of two lower than the num- 



(2010) due to differences in the calculation of the ionizing 
photon rate. For example, Dayal et al.| p009) uses the IMF 
mass range of [1,100] M rather than [0.1,100] M employed 
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Figure 3. The ionized region around the most massive object 
in L05 shaped by the density structure of the IGM close to the 
source. The degree of ionization of the gas is shown by the color 
bar at the top. 



in our simulations. Thus, their values would be comparable 

tO OUrS for /esc,ion,ISM = 0.04 and /esc, Lya, ISM = 0.6. 



4.3 CRASHa Input and Output 

The inputs for a CRASHa run are the gridded density, ve- 



locity and temperature fields (Sec. 4.1) and a source dis- 



tribution with photon rate and spectrum (Sec. 4.2). In our 



reference simulations, we adopt uniform initial fields of ion- 
ization, xion = null/fin = 0, and temperature, T = 100 K. 



© 2012 RAS, MNRAS 000, [Tp2] 



Effect of IGM on LAE detections during EoR 7 



A discussion on the different choices for these values is given 
in Section [6.31 The duration of the RT simulation is set to 
2 x 10 8 years. The source locations and luminosity were cal- 
culated as explained in the previous Section. Other simula- 
tion parameters, N p = 10 6 , iV P) Lya = 10 4 and Af em ,i = 500, 
were determined using resolution tests. 

Throughout the simulation, CRASHa collects photon 
packets exiting the box, recording their frequency, time, po- 
sition and directional information, which can be used to 
quantify observed properties of the source such as SB maps 
and spectra. Each cell on a side of the simulation cube emits 
photon packets in different directions. Just like the photon 
counts, the true distribution of angular directions of the pho- 
tons exiting from the cells are also sampled by the photon 
packets. Each direction points towards a different observer 
thus making each cell (3D) act as a pixel (2D) in the plane 
perpendicular to the direction of the observer. Here we de- 
fine the direction of the observer as the one perpendicular to 
the plane of the side of the cube. By this definition we get six 
different observer directions for each source. Both SB maps 
and spectra can be calculated for each observer direction. 

The SB SB em of a pixel of size p s due to a source of 
luminosity L at a distance d (in Euclidian geometry) is given 
by: 



SB e 



And 2 (Ps/d) 2 A-np 2 s 



(7) 



The photons, after travelling through the expanding uni- 
verse, arrive at the observer at z = 0. In a cosmological 
context, the observed SB SB oos of a source of luminosity L 
at redshift z — z s observed in a pixel of proper length p s is 
( Tolman||1934| ): 

SB ohs = j—^-X— X9 ,_„ 9 /-, ■ „ \a 77~ e ™zM ( 8 ) 



(Ps/Da) 2 47T P 2(1 + 2s )4 (l + Zs 



where Dl is the luminosity distance of the source and Da 
is its angular diameter distance. Therefore, to calculate the 
observed SB in a pixel for a specific observer, we need the 
number of photons in each cell travelling in the direction 
of the observer. The sampling of the angular distribution of 
photons depends on the number of photon packets in each 
cell. But the average number of photon packets in each cell 
of our simulation is too low (~ 12) for a smooth sampling 
of the underlying angular directional distribution. Thus we 
calculate the angular directional distribution using all the 
photon packets on a side of the box, by assuming that all 
cells on that side sample the same angular directional dis- 
tribution with respect to the true north for each cell, which 
is defined as the radial direction from the source to the cell. 

Figure [2] shows the sketch of the method. The star rep- 
resents the cell with the stellar source. The dotted lines rep- 
resent the radial direction connecting each cell on the edge 
of the box to the source (i.e. true north for each cell). The 
solid lines with arrows represent the directions of photon 
packets exiting such cells. The dashed lines from each cell 
to the observer represent the observer direction, which for 
z s = 7.7 is perpendicular to the side. The packets which are 
not scattered by the IGM follow the radial direction from the 
source (i.e. true north), while the scattered photon packets 
sample random directions. To calculate the SB ODS we need 
to estimate the number of photons reaching the observer at 
z = 0, i.e. traveling along the dashed lines. To calculate the 



underlying angular direction distribution, we grid the direc- 
tions of all photon packets (short solid arrows) in each cell, 
correcting for the different true north direction of each cell 
(dotted arrow). The gridded distribution of angles with re- 
spect to their true north in all the cells on a side are added 
and normalized using the total number of photon packets ex- 
iting the side. This normalized gridded angular distribution 
is then imposed on the photon count in each cell to obtain 
the number of photons in each pixel going in the direction 
of the observer. 

The total energy of photons travelling in a given direc- 
tion within a pixel in each second, L p i xe i, is given by: 

L 



lpixei — 



47ld 2 



X P s 



(9) 



which is calculated by adding up the energy of all the pho- 
tons in the pixel in the direction < p s /d. d is the proper 
distance of the cell to the source (at the centre of the sim- 
ulation box). Lpi X ei is converted to the units of SB, i.e. 
rad -2 , using the formula: 
2 



erg s 1 cm 2 



SB ern L D 



el X 



d z 



(10) 



This is then converted to the SB ODS dividing by (1 + z) as 
in equation (fsl) . The SB ODS is then converted to flux, F, as: 



F — SB hs x 



D 2 ' 

A 



(ii) 



where p s /Da is the pixel size of the observer in radians and 
Da is the angular diameter distance of the source. Adding 
up the flux in all pixels, we can obtain the observed source 
luminosity L Q bs as: 



L obs = 47rL>£ x F, 



(12) 



where Dl is the luminosity distance of the source. 

An important technical caveat is that our method for 
computing the SB differs from the often used 'next-event es- 
timator' approach which has been adopted in other studies 



e.g. Tasitsiomi|2006| Laursen, Sommer-Larsen, &; Andersen 



2009 



Zheng ltal.|20101|Kollmeier et al.|2010||Barnes et al. 



2011 ). Idealized test cases though show that the method per- 



fectly reproduces the SB profiles of the Rybicki-Loeb halo 
(Loeb &; Rybicki 1999), and also those produced in cases 



where a spherical symmetric ionized bubble is carved out 
of the neutral expanding IGM (as in Dijkstra & Wyithe 
2010). This approach though is less appropriate for cases 
where scattering occurs in only a few (neutral) clumps em- 
bedded in an otherwise fully transparent medium. However, 
the method works well in test cases similar to the config- 
urations we are studying in our simulations, and thus we 
consider the main results of this study to be robust. 



5 RESULTS 

Using the pipeline outlined in the previous Section, we sim- 
ulate ~ 130 LAEs from L05-30 to evenly sample the DM 
mass range 10 9-11 M©. Table [2] contains the properties of 
the haloes used in this work. The simulations were run with 
our default parameters (see Sec. |4.3| ). Dependence of the 
results on the choice of the parameters is discussed in Sec- 
tion [6] Before trying to understand the properties and sta- 
tistical trends of 130 galaxies, we focus on understanding 
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Figure 4. Contour plot showing the six different SB maps for the same object as obtained from the six sides of the simulation cube. 
The object in this simulation is the most massive object in L05, the same one as in Figure [T] The SB values in the pixels are shown in 
the color bar at the top. See text for details. 




arc sec 



Figure 5. A cross section of the SB profile for the maps in the 
top row of Figure ^] The lines are for SB values along the y axis 
(i.e. x values are 0) of the maps. The solid red, dashed green and 
dotted blue lines refer to the top- left, top-middle and top-right 
maps, respectively. 



how the IGM close to the objects affects the appearance of 
a single LAE. We investigate the ionization structure in the 
IGM and how it affects the SB maps. Having understood 
the effects on an individual galaxy, we study the behaviour 
of the whole sample and the trends shown in observability 
and Lya escape fractions due to the IGM. 



5.1 Behaviour of an Individual Galaxy 

To understand how the density field in the IGM around the 
source affects the shape of the ionized region, which in turn 
will reflect on the propagation of the Lya photons, we look 
at the ionization structure around the most massive dark 
matter halo in L05 (same galaxy as in Fig. [I) of DM mass 
of 2.6 x 10 10 M and stellar mass of 3.6 x 10 8 M . Fig- 
ure [3] shows the mid-plane of the ionization structure at the 
end of the RT simulation. We can see that the ionized re- 
gion is not spherical and its edges are shaped by the high 
density regions (e.g. filaments) in the IGM. Lya photons 
propagating through the ionized gas will travel a large dis- 
tance unscattered while being redshifted away from the Lya 
rest wavelength, thus improving the chances of escape from 
the neutral IGM beyond the ionized region. Lya photons 
encountering a high density filament instead will undergo 
several scatterings before being able to eventually escape 
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Figure 6. Stellar mass of the objects plotted against their DM 
mass. The objects from different simulations are marked as - L05 
(red crosses), L10 (green squares), L20 (blue circles) and L30 
(pink triangles). The best fit line (dashed) is obtained using only 
haloes from L05 and L10, and it has a log- log slope of 1.12 with 
a standard deviation cr=0.13 dex. The a region is shaded in light 
grey. 
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Figure 7. Observed Lyo: luminosity of the objects plotted against 
their DM mass. The objects from different simulations are marked 
as - L05 (red crosses), L10 (green squares), L20 (blue circles) and 
L30 (pink triangles). The best fit line (dashed) is obtained using 
only haloes from L05 and L10, and it has a log- log slope of ~ 1.12 
with a standard deviation cr=0.17 dex. The a region is shaded in 
light grey. 



the simulation box. Because of the different paths followed 
by the photons, we expect different SB distributions and 
spectra, depending on the viewing direction. This can be 
clearly seen in Figure [4] where the SB maps obtained from 
the six sides of the simulation box are shown for the same 
object. 

The differences in the maps reflect the structure in the 
IGM. In general, the lesser the photons are scattered before 
reaching the edge of the box, the more concentrate the SB 
profile is. In this case, the central pixels also have a much 
higher value of SB compared to cases in which the photons 
are scattered by high density neutral gas. Note that, as a 
result of the very inhomogeneous structure of the IGM, the 
maximum SB value in the image differs by more than an 
order of magnitude for different lines-of-sight. This is more 
evident in Figure [5] which shows a cross section of the SB 
maps for the three directions shown in the top row of Fig- 
ure [4] Plotted are the SB values in the map for an x axis 
value of against the y axis of the map. As noted earlier, 
the difference in the value of the central pixels in the maps 
is more than an order of magnitude for different viewing 
directions of the same object. 

The detect ability of the objects simulated in this work 
depends on the SB thresholds SBth set by different obser- 
vational campaigns. Since the lowest SB value of a pixel in 
the maps is SB m i n — 10 -23 erg s _1 cm -2 arcsec -2 , only for 
such a small (or lower) value of S'.Bth would it be possible 
to detect all the flux from this object. At this SBth level 
the observed luminosities, calculated using Equations [TT] 
and [12] for the source shown in the previous Figures are 
in the range L Q b s = 4 — 12 x 10 41 erg s _1 . For reference, 
the input luminosity is L^ a = 9.2 x 10 41 erg s _1 . This 
would lead to an escape fraction / eS c,Lya,iGM = L Q h s /L^y a 
in the range 0.43 — 1.3 (see below for a discussion). If in- 
stead e.g. SB t h — 3 x 10 -21 erg s _1 cm -2 arcsec -2 , only 
two maps would result in a detection, with luminosities 
L obs (> SBth) = 8xl0 40 erg s" 1 and 7.9 xlO 39 erg s" 1 , com- 
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Figure 8. Lyo; photon escape fractions due to the IGM, 
/esc,Lya,lGM 5 plotted against DM mass for objects from L05 (red 
crosses), L10 (green squares), L20 (blue circles) and L30 (pink 
triangles) simulations. The best fit line has a log-normal slope of 
-0.01, with mean / eS c,LyQ,lGM of ~ 0.73 and a standard deviation 
cr=0.18. The a region is shaded in light grey. 

pared to L obs (> SB min ) = 1.2 x 10 42 erg s" 1 and 8 x 10 41 
erg s _1 , respectively. Note that if observations had a thresh- 
old SB at this level only a few percent of the photons would 
be observed leading to a very low effective Lya escape frac- 
tion / e f c , Ly « ;IGM = L obs (>_SB th )/L?° = 0.007 - 0.08. For 
SBth = 10 21 erg s 1 cm 2 arcsec 2 , detections would be 
made in four maps with an observed luminosity in the range 
L bs(> SBth) = 2 x 10 39 - 5 x 10 41 erg s~\ corresponding 
to /efcLy^iGM = 0.002 - 0.6 for Lff y c a = 9.2 x 10 41 erg b" 1 . 

Some lines-of-sight can have an effective escape fraction 
> 1. These are lines-of-sight through voids, which in addition 
to the isotropic flux of photons from the central source, also 
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Figure 9. Observed Lyoi luminosity of the objects from all simulations plotted against their DM mass. Symbols refer to different SB 
cuts, SB^: 10 -25 erg s _1 cm -2 arcsec -2 (red crosses), 10 -21 erg s _1 cm -2 arcsec -2 (green squares), 10 -20 erg s -1 cm -2 arcsec -2 
(blue circles) and 5 x 10 — 20 erg s — 1 cm -2 arcsec -2 (pink triangle). 



get an additional contribution of photons which have been 
scattered off from other lines-of-sight. A Lya photon encoun- 
tering the surface of a high density region (like a filament 
close to the source) has a high probability of scattering to- 
wards the lower density gas around it. Once in a line-of-sight 
through a void, the photon travels longer distances with- 
out scattering while redshifting out of resonance, improving 
the probability of escape towards the observer. Thus lines- 
of-sight through voids are preferred over the ones through 
high density regions and get Lya photon contributions from 
other lines-of-sight, leading to a more than average photon 
escape and an occasional effective escape fraction > 1. Note 
that while in the case of absorption photons removed from 
a line-of-sight do not contribute to any other line-of-sight 
leading to an escape fraction always < 1, in the case of scat- 
tering photons removed from a line-of-sight appear in an- 
other one leading to escape fractions which can be also > 1. 
This makes a 3D treatment of the Lya photon scattering of 
particular relevance for a proper LAE modelling. 

5.2 Statistical Trends 

Next we discuss how the distribution in SB affects the sta- 
tistical properties of a large sample of LAEs. The difference 
in SB profiles in fact leads also to a spread in the Lya lumi- 
nosity of the same object observed from different directions, 
as seen earlier. The objects in our sample have an intrinsic 
correlation between the stellar and the DM mass, with a log- 
log slope of 1.12 and a standard deviation cr=0.13 dex (see 
Fig. |6j). The best fit line (dashed) has been obtained using 
only haloes from L05 and L10. In fact, due to low particle 
resolution in the larger boxes, some of the haloes selected 
in L20 (blue circles) and L30 (pink triangles) have a stel- 



lar mass which lies below the expected trend. A correlation 
similar to the one shown in Figure [5] is also expected to exist 
between the input Lya luminosity of the source, L^fya? an d 
the DM mass. Any additional scatter in the observed Lya 
luminosity L Q bs is instead due to the effect of the IGM. 

Figure [7] shows L bs of the full sample calculated for 
each of the six sides of the box plotted against the DM mass 
for each of the simulated LAEs. As in Figure [6] there is a 
strong correlation with a best fit log- log slope of 1.12, which 
has been calculated using only haloes from L05 and L10. 
The scatter now though is 0.17 dex, larger than the intrinsic 
one, and it is induced by the structure in the IGM. 

An alternative way to look at the scatter due to the 
IGM is to plot the escape fraction / eS c,Lya,iGM against the 
DM mass of the objects (Fig. [8j). This removes the scatter 
present in the intrinsic luminosity L^y^. As we can see, there 
is no strong correlation between / eS c,L y a,iGM and the DM 
mass because the scatter due to the IGM structure dom- 
inates. The average escape fraction has a nearly constant 
value of rsj 0.73 with a a scatter of 0.18. Some lines-of-sight 
have escape fractions greater than 1, as explained earlier. 
But the statistical mean of the six observed IGM escape 
fractions for all the objects in the sample is 0.73. We would 
get an average escape fraction of 1 only if we observed all the 
directions for the same object and calculated the mean of all 
the lines-of-sight. The average of the six viewing directions 
is < 1, showing that scattering removes a large fraction of 
the flux for most lines-of-sight. Most of the flux from the ob- 
ject escapes along a small fraction of lines-of-sight through 
large voids, showing again the importance of IGM struc- 
ture close to the source for the observability of LAEs. The 
spread in the correlation decreases for higher DM masses. 
This could be due to two reasons - a more uniform environ- 
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Figure 11. Surface brightness maps for one of the sides of the cube for runs with different /esc, ion, ISM- Plotted are the maps from the 
simulations with /esc, ion, ISM = 0.02 (left; reference case), 0.5 (middle) and 0.99 (right). See text for details. 
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Figure 10. Spectrum of all the Lya photons exiting the sim- 
ulation box for runs with different /esc, ion, ISM- The lines refer 
to /esc, ion .ISM — 0.02 (red solid line; reference case), 0.5 (green 
dashed) and 0.99 (blue dotted). The vertical black line marks the 
Lyo; rest-frame wavelength of 1215.67 A. 
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Figure 12. Spectrum of all the Lya photons exiting the simu- 
lation box for runs with different /esc, Lye*, ISM- The lines refer to 
/esc,Lya,lSM = 0.01 (green dashed line), 0.3 (red solid line; refer- 
ence case) and 0.99 (blue dotted). The vertical black line marks 
the Lyo; rest-frame wavelength of 1215.67 A. 



ment for haloes of higher masses leading to lesser scatter; or 
the resolution effects caused by the larger physical size of the 
grid cells (oc r2oo) for higher mass haloes compared to low 
mass ones. But the general consistency of / e sc,L y «,iGM values 
in L20 and L30 with L05 seems to indicate that resolution 
issues might be less important. Higher resolution RT sim- 
ulations of larger comoving volumes are needed to confirm 
this effect, which is beyond the scope of this study. 

Next we investigate the detectability of this dataset by 
plotting in Figure [9] the observed luminosity of the objects 
calculated from the six sides of the simulation boxes for dif- 
ferent SB thresholds, SBth- Only one object is detected for 
SB t h ^ 5 x 10 -20 erg s -1 cm -2 arcsec -2 (pink triangle). 
Note that this is not the most massive object neither in the 
full sample nor in the simulation box L05. Nevertheless, it 
looks as the brightest due to the density and velocity fields 
in the surrounding IGM, which lead to a very concentrated 



SB profile and ease the detectability of the source. From this 
example it is clear that the structure in the IGM plays an im- 
portant role in the detectability of objects by shaping their 
SB profile and only a full 3D RT approach can properly cap- 
ture its effect. As we lower the SB threshold, more objects 
are detected and a large fraction of the luminosity is ac- 
counted for. By SBth = 10 -20 erg s -1 cm -2 arcsec -2 (blue 
circles), seven objects are detected covering a mass range of 
two orders of magnitude. Even though there is a correlation 
of luminosity with DM mass apparent in the data, a large 
scatter is present at small masses. Decreasing the threshold 
by another dex, leads to the detection of almost all objects 
with a scatter in luminosity spanning three dex in each mass 
range. A SB level of SBth = 10 -25 erg s -1 cm -2 arcsec -2 
is needed to observe the total flux from all the sources (the 
same as Fig. FFb. The SB values calculated here depend on 
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the choice of a number of different parameters, which will 
be further discussed in the next Section. 

We caution the reader that current observational SB 
thresholds are significantly higher than the 'typical' SB lev- 
els of Lya radiation scattered in the IGM as found in our 
simulations. This leads to an effective escape fraction due to 
the IGM (also see Laursen et al. 2011) which is very low (a 
few %). In addition, the observed SB is SB Q h s oc (1 + z)~ 4 
(see Eq. [8|, making deep detections increasingly difficult at 
higher redshifts. Current observational thresholds for nar- 
row band surveys are 2.64 x 10 -18 erg s _1 cm -2 arcsec -2 
at z = 5.7 in SXDS (Subaru/XM M- Newton Deep S u rvey; 
M. Ouchi; privat e communication; [Zheng et al.||2010 ). Stei- 
del et al.| ( |201l| ) achieved deep SB limits of 2.5 x lCT iy 
erg s~ cm - arcsec -2 for Lyman Break Galaxies (LBGs) 
at z rsj 3 which was obtained through stacking of 92 im- 
ages. Their total observed Lya flux translates to an escape 
fraction of Lya of ~ 0.25 — 0.5 (see Dijkstra & Hultman 
Kramer|[2012t . This is in agreement with Figure |9j Ob- 
taining these SB thresholds at higher redshifts is a very 
difficult task. JWST, with its predictecQ SB th ~ 10~ 19 
erg s _1 cm -2 arcsec -2 , would be an important step for the 
deep detections of LAEs. 



6 PARAMETER STUDY 

In this Section we study the dependence of our results on the 
different parameters adopted. To this aim we use the biggest 
object in L05 with a dark matter mass of 2.5 x 10 10 M©. 
The reference values are those described in Section T4.3I 

A minor contribution to the observed flux comes from 
the photons emitted by the recombining gas in the IGM, 
which has been estimated to be ~ 3%. This is not signif- 
icant compared to the flux from the central source, but it 
should be kept in mind that this flux is not affected by the 
absorption in the ISM and thus has a higher escape prob- 
ability. Since it is not concentrated as a point source, its 
contribution to the detect ability of the object is negligible; 
rather, it might reduce the strength of the point source by 
spreading the flux to the scattered component of the angular 
directional distribution. We do not include the re-emission 
photons in our parameter study. 



6.1 Escape Fraction 

As already discussed, we need to specify the value of both 
ionizing, / eS c,ion,iSM, and Lya, / e sc,L yc *,iSM, escape fractions 
from the ISM. The amount of ionizing photons reaching the 
IGM is controlled by / eS c,ion,iSM, thus determining the ex- 
tent of the ionized region through which the Lya photons 
propagate. The value of the escape fraction also affects the 
production of Lya photons by recombinations in the ISM. In 
our reference simulations, the box around the object has a 
size of ~ 35 r2oo comoving. Because this is not large enough 
to contain the HII region produced by /esc, ion, ism = 0.99, 
for this parameter study we use an object with DM mass of 



3 From http://www.stsci.edu/jwst/science/sensitivity for 3cr de- 
tection in 10 5 s with NIRSpec IFU mode at 0.1 x 3 arcsec 2 reso- 
lution. 
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Figure 13. Spectrum of all the Lya photons exiting the simula- 
tion box for different velocity shift of the input Lyo; monochro- 
matic spectrum. The lines refer to a velocity shift of (red solid; 
reference case), 200 (green dashed), 400 (blue short dashed) and 
800 (pink dotted) km s _1 . 



2.29 x 10 10 M from L10. We extract a cube with a side of 
~ 117 r2oo comoving, which is gridded to a 256 3 grid. 

To quantify the effect of changing the value of 
/esc, ion, ism on the spectrum, in Figure [l0| we show the Lya 
spectrum of the photons exiting the box for simulations with 
/esc, ion, ism = 0.02,0.5,0.99, while /esc, Lya, ism is set to the 
default value of 0.3. Note that the photon packets are al- 
ways blueshifted back to the source position from the edge 
of the box. Also note that the spectrum is for all the pho- 
tons exiting the box throughout the simulation time. The 
Lya photon count increases with decreasing /esc, ion, ism (see 
Eq.[5]), as recombination in the ISM is the dominant contri- 
bution to Lya photon production. It also affects the shape 
of the spectrum by moving its peak closer to 1215.67 A, as 
the photons escape with less scatter. This can also be seen 
from the SB maps in Figure [TT] The higher the /esc, ion, ism 
is, the more concentrated is the SB profile. But due to the 
lower production of Lya photons in the ISM, the SB values 
in each pixel decrease for higher /esc, ion, ism - Thus, increasing 
/esc, ion, ism reduces the overall detect ability of the source at 
a specific SBth , but improves its detectability at lower SBth 
with respect to a dimmer galaxy with lower / esc , ion, ISM giv- 
ing similar total Lya fluxes. 

Direct estimation of /esc, ion, ism from observations is a 



very difficult task (e.g. 


Bland-Hawthorn & Putman 


2001 


|Steidel, Pettini, & Ade 


lberger||2001| |Shapley et al. 


2006 



into modelling the correlation between halo properties and 
ionizing escape fraction from the galaxy (e.g. Dove, Shull, 



Ferrara||2000| |Wood Loeb 20 00] |Ricotti fe Shull||200^ 
|Ciardi, Bianchi, fe Ferrara||2002| |Clarke fe Oey||2002| |Fu- 
jita et al. 2003; Razou mov fc Sommer-Larsen|2006| |Gn edin, 
Kravtsov, fe Chen||2008| |Wise fc Cenl|2009| |Yajima et al. 
2009 Yajima, Choi, &; Nagamine|201l[|Paardekooper et al. 



2011 Yajima et al.||20lT ). Estimated escape fraction values 
range from 10~ to 1 depending, among others, on the mass, 
redshift, dust and gas distribution and star formation rate, 
with different studies giving different values. Since the cor- 
relation between halo properties and /esc, ion, ism is not well 
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Figure 14. Surface brightness maps for one of the sides of the simulation box and for different velocity shift of the input Lya monochro- 
matic spectrum, i.e. (bottom left panel), 200 (bottom right), 400 (top left) and 800 (top right) km s _1 . 




Figure 15. Cross section of the SB profile for the maps in Fig- 
ure^] Plotted are the SB values along the y axis for an x axis 
value of 0. The lines refer to a velocity shift of the input Lycn 
monochromatic spectrum of km s — 1 (red solid; reference case), 
200 km s _1 (green dashed), 400 km s _1 (blue short dashed) and 
800 km s" 1 (pink dotted). 



known, we choose our default value to be f ei 
Gnedin, Kravtsov, Chen| ( pOO^ 



q,ism — 0.02, 



of Lya photons due to the ISM, / eS c,Lya,iSM, which controls 
the fraction of Lya photons reaching the IGM, as detailed 
in Equation [6] We use the same object and box of the pre- 
vious tests, with /esc,L y a,iSM = 0.01,0.3 and 0.99. Figure [l2| 
shows the spectrum of the photons exiting the cube, whose 
intensity increases with increasing /esc,L y( *,iSM as more pho- 
tons reach the IGM. On the other hand, the shape of the 
spectrum is not altered. Similarly, the morphology of the SB 
profile is not affected. 



The escape fraction of Lya photons from the ISM has 



been estimated in a large number of studies (e.g. Le Del 



liou et al.|2005 |2006| 


Dave, Finlator, & Oppenheimer 


2006 


Nagamine et al.||2010 


|Ouchi et al. ||2008 |Dayal et al. 


2009 


Kobayashi, Totani, & Nagashima 2007; Laursen, Sommer- 


Larsen, & Andersen|2009 ; Steidel et al.|2011||Schaerer et al. 


2011||Forero-Romero et al.|2011| |Yajima et al.|2011). 


Since 



the amount and type of dust in the ISM of a galaxy and its 
correlation with other halo properties is highly unknown, we 
choose a reference value of 0.3. This adopted value is sim- 



We have then tested the effect of the escape fraction 



ilar to the value that was required by Dayal et al. (2009) 
in order to reproduce the observed Lya luminosity func- 
tions at z = 5.7 (note that this requirement assumes that 
the IGM transmits ~ 50% of all photons: for a more trans- 
parent/opaque IGM, the required escape fraction would be 
lower /higher). 
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Figure 17. Spectrum of all the Lya photons exiting the simula- 
tion box for runs with initial volume averaged ionization fraction 
Xion=0 (red solid line; reference case), 0.5 (green dashed) and 0.89 
(blue dotted). 
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Figure 19. Same as Figure^] but including scattering of Lyo; 
photons by the column of neutral gas beyond the edge of the 
cube. The y-axis has the same scale as Figure |17| The lines are 
for a volume averaged ionization fraction Xi on =0 (red solid line; 
reference case), 0.5 (green dashed) and 0.89 (blue dotted). 



6.2 Effect of Input Lya Spectrum 

The RT of Lya photons through the IGM is affected by their 
wavelength as they escape the ISM. In our default runs, we 
use a monochromatic line at 1215.67 A, but we can expect 
the spectral shape of the photons emitted by the stellar 
sources to be distorted by their interaction with the ISM 
before entering the IGM. In fact, observations of LAEs at 
z rsj 2 — 3 show a non-monochromatic spectru m (|Verhamme 
et al.||2008). Other authors (e.g. payed et aL||2009{ |Zheng 
Dijk stra &; Wyithe||2010 ) use gaussian profiles 



et al. 2010 



with a width given by a fraction of the circular velocity at 
the virial radius of the DM halo, with wider profiles leading 
to an easier escape of the Lya photons from the surround- 
ing medium. A lot of effort has gone into modelling more 
accurately the Lya line profile for a range of parameters 
(Sch aerer et al.||2011| , but the exact shape of the spectrum 
and its connection to the galaxy properties are still unclear. 
Because of all these uncertainties, we choose to investigate 
the effect of the input Lya spectrum simply by shifting a 
monochromatic spectrum. 



Figure [13] shows the spectrum of all the Lya photons ex- 
iting the box for different values of the input Lya monochro- 
matic photons, i.e. with velocity shifts of (reference case), 
200, 400 and 800 km s _1 . Since the scattering cross sec- 
tion is oc Az/ -2 , the larger the shift, the lesser the scat- 
ter suffered by the photons. The velocity field in the IGM 
also plays a role in determining the scattering probability 
of these photons. Since the velocity distribution of the IGM 
around the object has a complex structure, it is not possible 
to uniquely predict the velocity shift at which the photons 
stop being scattered. For our test object, as the photons are 
shifted by 200 km s _1 , the scattering is reduced by a factor 
of two, which in turn is reflected in a reduction of photons 
at 1215.67 A. For a shift of 400 km s~\ the scattering is 
reduced by an order of magnitude. Monochromatic spectra 
shifted by > 800 km s _1 get hardly scattered. This also af- 
fects the SB profiles, as shown in Figure [l4| where the SB 
in the four cases is plotted for one of the sides of the box 
(top-left map in Figure [4]). As the velocity shift increases, 
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Figure 18. Surface brightness maps for one of the sides of the box for runs with initial volume averaged ionization fraction Xi on = (left 
panel; reference case), 0.5 (middle) and 0.89 (right). 



the probability of scatter of photons decreases, leading to a 
more concentrated SB profile. In addition, only a small frac- 
tion of the photons are scattered, forming a low SB halo. 
This can be seen more clearly in Figure [15] where the cross 
section of the maps is shown. The larger the velocity shift 
is, the higher is the SB value in the central pixel, improv- 
ing the detect ability of the source. Also, a larger fraction of 
the flux escapes from this pixel. For example, 0.004, 0.1, 1.5 
and 4 per cent of the total flux escapes through this pixel 
for the cases of no velocity shift, a velocity shift of 200, 400 
and 800 km s _1 , respectively. Therefore, the initial spectra 
of Lya photons (i.e. the shape induced by the processing of 
the Lya photons in the ISM) has a very important effect in 
determining the SB profile of the object and on the observ- 
ability of LAEs. Quantifying this in more details is beyond 
the scope of this study and will be investigated elsewhere. 

From Lya spectrum models of Sch aerer et al.| ( |2011| ), 
the peak of the spectra of photons escaping from the galaxy 
is at wavelengths where the velocity shift is equivalent to 
~ 1.5 — 2 times the outflow velocities. The typical outflow 
velocit ies estimated from observations by Verham me et al7| 
(2008) is 150 — 200 km s _1 , thus showing a Lya spec- 
trum peak at 200 — 400 km s _1 . With a simplifying assump- 
tion that almost all the flux enters the IGM at these wave- 
length, we can estimate the expected changes in observed 
luminosity at different SB t h due to this velocity shift. At 
^S^th = 10 -21 erg s _1 cm -2 arcsec -2 , for a velocity shift 



of 200 (400) km s" 1 , / e f c , Lya)] 



0.01 (0.1), compared to 



7 x 10~ in our reference case with no velocity shift. Thus 
detectability of a LAE at a specific SBth can improve by 
orders of magnitude for realistic velocit y shifts of 200 — 400 
km s _1 even with Xi on = (see also |Dijkstra &; Wyithe| 
2010). Therefore, understanding the effect of the ISM of the 



galaxies along with their outflows is a crucial input for better 
modelling of LAEs. This topic will be addressed elsewhere 
in further detail. 



4 Note that a specific SB threshold might select more than 1 
pixel. 



6.3 Effect of IGM Ionization 

Another important factor determining the observability of 
LAEs is the ionization level of the IGM outside the fully 
ionized region created by the source. The initial ionization 
fraction of our simulations is set to a uniform default value of 
xion = for convenience, but simulations of the reionization 
process indicate that most of the IGM is substantially ion- 



ized at the redshifts studied in this paper (e.g. Ciardi et al 
|2Q11| ). The IGM ionization level is determined by a combi- 
nation of a background radiation and local sources. Here we 
study how these two contributions change the observability 
of our test source. 



6.3.1 Uniform Ionization Distribution 

To understand the effect of the IGM ionization level on the 
Lya RT, we repeat the reference simulation with an initial 
volume averaged ionization fraction of Xi on = 0.5 and 0.89. 
The first value corresponds to an epoch when only half of 
the IGM is ionized, while the second is the value obtained 
from the reionization model £1.2-al.8 described in Ciardi et 
ED (12011), which is designed to be consistent with existing 
observational constraints such as the Thomson scattering 
optical depth and the photo- ionization rate at z < 6. 

Figure [l6| shows a slice cut through the ionization struc- 
ture at the end of the simulation for the three cases described 
above. As expected, the higher Xi on is, the larger is the fully 
ionized region produced by the source, although the differ- 
ences are not dramatic because the high density filaments 
surrounding the source effectively confine the fully ionized 
region in all cases and determine its shape. Note that while 
the edges of the ionized region are controlled by the high 
density filaments in the IGM, the inclusion of ionizing flux 
from galaxies in those high density clumps could lead to a 



different ionization structure (as shown in Sec. 6.3.2). 

As a larger ionized region results in less scatter for the 
Lya photons, we expect an easier escape in this case. This 
is shown in Figure |17[ where the spectrum of all the Lya 
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Figure 20. A slice through the ionization structure for RT simulations showing the effect of clustering of sources. The ionized region 
is produced by the central source alone (left panel), all the sources in the simulation box (case C, middle panel) and the central source 
alone in a uniform partially ionized IGM (case 14, right panel). The final volume averaged ionization fraction in the middle and right 
panel is the same. See text for further details. 
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Figure 21. Surface brightness maps for one of the sides of the simulation box for the same simulations of Figure [20] 



photons exiting the simulation box is plotted. Due to the 
lower number of scatterings undergone by the Lya photons 
for higher values of Xi on , the peak shifts to bluer wave- 
lengths and the photon count at the peak increases while 
reducing the width of the line profile. The SB maps for 
one of the sides of the box for the different Xi on values is 
shown in Figure[l8] As Xi on increases, the SB profile becomes 
more concentrated while improving detectability. From the 
maps, we can infer that observations at SB t h = 3 x 10 -21 
erg s _1 cm -2 arcsec -2 would appear as a point source to the 
observer. For this viewing direction, the luminosity (effective 
escape fraction) of the point source detected for Xi on = 0.5 
and 0.89 is 8.4 x 10 38 erg s" 1 (0.001) and 3.7 x 10 40 erg s" 1 
(0.04), respectively, compared to a non detection at Xi on = 0. 
More specifically, if we were to consider the maps from all 
six faces of the cube, at Xi on = 0.89 the object would be 
detected as a point source from all the six directions with 
/esc,Lya,iGM m the range 0.004 — 0.6, compared to only two 



detections as point source in six maps at Xi on = with 
/es^LyajGM = 0.008 and 0.09. Thus at high ionization frac- 
tions most of the objects can be detected as point sources 
at lower SB thresholds compared to our reference case. Our 



results appear similar to Zheng et al. 



(2010) who find that 



only a small fraction 



even at a mean neutral fraction of 10" 
of the photons (8-33 % of the intrinsic luminosity) escape as 
a point source. 



Because this paper is primarily focused on the effect 
of the IGM in the immediate surroundings of a source of 
Lya photons, we do not follow the propagation of the pho- 
tons to the observer. To do so, we would need much larger 
boxes, losing the resolution necessary to properly account 
for the scales we are mainly interested in here. Nevertheless, 
the IGM outside our simulation box is bound to scatter the 
photons in an amount dependent on their frequency when 
they escape the box. The optical depth r due to scattering 
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away from the line- of- sight to the observer (at z — 0) for a 
LAE at z — z s is given in Equations ^ and [2] 

For xhi = 1 at z s = 7.7, the optical depth r ~ 1 for 
~ 1400 km s _1 , corresponding to a comoving distance 
of 8.6 /i _1 Mpc from the source which is much larger than 
our cube sizes. Beyond this distance from the source the Lya 
photons have redshifted far enough from resonance such that 
the neutral IGM does not further affect the Lya SB profiles. 
However, in some cases analyzed in this work - in particular 
low mass objects for which smaller boxes are extracted - 
the photons emitted by the source do not redshift far from 
the line resonance when they reach the edge of the box. The 
neutral gas outside the box can then further modify the 
predicted Lya SB profile. 

As a simple estimate of the scattering in neutral gas 
beyond the boundaries of the cube we can use Equation [2] 
The effect of this scattering on the spectra of photons es- 
caping the box for our test object is shown in Figure |19| 
For Xion = 0, the r integrated over the whole spectrum is 
2.9, but as the ionization increases, the r decreases to 1.8 at 
Xion = 0.5 and to 0.55 at Xi on = 0.89, leading to more pho- 
tons being observed. The scattering beyond the cube can 
therefore be significant for small Xi on , resulting in lower SB 
values per pixel extending over a larger angular size. Even 
deeper SB cuts would then be required to observe the to- 
tal luminosity of the object. In contrast, for high Xi on scat- 
tering beyond the cube can be safely ignored, and the SB 
profile would retain the point source like appearance seen 
in Figure [T8] The qualitative results described in the previ- 
ous sections (e.g. a spread and increment in /escfLya,iGM f° r 
lower SBth) would still hold. Another factor which reduces 
scattering outside t he s imulation box is the velocity shift 
discussed in Section 6.2 A shift of 200 (400) km s _1 reduces 
the effective r to 2.8 (2.3) even for Xi on = 0. This implies 
that the predicted steepening of the SB profiles with Xi on 
(see Sec. [7] below) and velocity off-set are (conservatively) 
underestimated by our calculations. 

We point out that constraints on the reionization his- 
tory that are provided by current observations of the Lya 
forest and the cosmic microwave background, as well as 
theoretical investigations, suggest that the ionization frac- 
tion at z = 7.7 was mostly complete (i.e. well in excess of 
0.50; see e.g. F ig. 11 of |Pritchard et al.|2010| or Fig. 4 



of Ciardi et al. 2011). This suggest that for realistic con- 
straints on Xi on at z — 7.7 the scatter beyond the cube is 
not significant, and so we do not include this extra scatter 
in the rest of our analysis. 

6. 3. 2 Clustering of Nearby Sources 

While a roughly uniform ionization fraction could be in- 
duced by an ionizing background, photons from neighboring 
haloes clustered around a source would induce further fluc- 
tuations in the gas ionization structure. To investigate the 
impact this has on our reference simulation, we propagate 
ionizing photons from all the stellar sources present in the 
cube (referred to as 'case C for clustering), rather than only 
from the central source. Figure [20] shows a slice through 
the ionization structure obtained in case C (middle panel) 
along with that from the reference simulation (left panel). As 
expected, the ionization bubble around the object is much 
larger and complex when the effect of all the sources is con- 



sidered. In case C, the final value for the volume averaged 
ionization fraction is Xi on = 0.340, compared to 0.047 of the 
reference case. 

To understand how the ionization structure due to clus- 
tered sources affects the propagation of Lya photons com- 
pared to a uniform ionization, we do another test run (re- 
ferred to as 'case W for uniform) producing the same final 
volume averaged ionization fraction of case C, i.e. 0.340, but 
this is induced by photons emitted only from the central 
source embedded in a gas with a uniform initial ionization 
fraction of 0.293. Case IA is shown in the right panel of Fig- 
ure[20]and exhibits a very different ionization structure com- 
pared to case C. While the region with fully ionized gas is 
smaller in case U, the situation is reversed in the region be- 
yond the bubble with some fully neutral gas still present in 
case C. Note that in all cases the Lya photons are emitted 
only by the central source. 

Figure shows the SB maps for the three cases dis- 
cussed above (for the top right panel in Fig. [4]). Both case 
C and case U show a more concentrated SB profile com- 
pared to the reference run as expected for a lower Xi on , 
but the details in the SB maps differ due to the distribu- 
tion of neutral gas. In particular, all the maps (from the 
six faces of the cube) of case IA have a more concentrated 
structure compared to case C due to the higher ionization 
level of the gas outside the fully ionized bubble. In fact, 
since the ionization bubble is not large enough for the Lya 
photons to be redshifted out of resonance when they reach 
its border, the neutral gas distribution beyond it is particu- 
larly important for the determination of the corresponding 
SB maps. Because the concentration in the SB improves 
detectability, case IA is expected to favor observations of 
LAEs. At SBth = 10 -21 erg s _1 cm -2 arcsec -2 , case C has 
/esc, Lya, igm °° 0.001, while case IA has ~ 0.004, compared to 
the reference case of ~ 0.0007. 

It might seem that clustering of neighbouring sources 
does not substantially change the observability of the LAE, 
but because the flux is redistributed in different direc- 
tions due to the different ionization structure in the IGM, 
/esc, Lya, igm can vary by orders of magnitude depending on 
the observing direction. For the same SB t h, another direc- 
tion (top right panel in Fig. Ul has /es^, Lya, igm = 1-5, 0.7 
and 0.6 for case C, case IA and the reference case, respectively. 
For reference, the input luminosity is Lf^ a — 9.2 x 10 41 
erg s _1 . Thus it should be kept in mind that even though 
both neighbouring sources and uniform background flux im- 
proves detectability of LAEs by increasing the mean ioniza- 
tion fraction of IGM in the region, the detailed ionization 
structure around a source plays a crucial role in determining 
the SB profiles in different directions. 



7 EFFECTS ON ESTIMATES OF THE 
REIONIZATION HISTORY 

Keeping in mind the discussions in the previous Sections, 
we now turn to study how structure in the IGM affects the 
methods which use the SB profiles of LAEs to estimate the 
mean ionization fraction of the universe. One of the most 
popular methods is to calculate the luminosity function of 



LAEs and determine its redshift evolution (e.g. Malhotra & 



Rhoads|2004| |Dayal et al.|2009| [Kashikawa et al.|2011[ ). Us- 
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Figure 22. Luminosity functions calculated for the 45 most massive objects in L05 for different SB cuts. The runs are for initial 
uniform volume averaged ionization fraction Xi on =0 (thick red), Xi on = 0.5 (medium thick green) and Xi on = 0.89 (thin blue). The SB 
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ing analytic models to describe the reionization history, the 
above authors estimate the evolution of r (employing equa- 
tions similar to Eq. [5]). The intrinsic LF of LAEs is then 
calculated from halo mass functions, assuming that the lu- 
minosity of an object is proportional to its mass and that 
it is attenuated by e _r due IGM absorption/scattering. A 
comparison with the observed LF provides a constraint on 
the actual r and thus on the gas neutral fraction in the IGM. 
But these methods rely on the assumption that observations 
of LAEs obtain the total flux from the source detected as a 
point source. Some recent observations though have detected 
faint extended SB profiles ( |Steidel et al.|2 011), as produced 
in our simulations. As we have seen, the probability of de- 
tection of a LAE strongly depends on SB t h-> especially for 
low xion. Here we investigate how a different choice of SB t h 
affects the LAE LFs for different values of Xi on , aiming at a 
qualitative assessment, while a more detailed and quantita- 
tive comparison to observations is part of a future study. 

We model the 45 most massive objects in L05 using 
Xion=0, 0.5 and 0.89 as described in Section |6.3.1| For 
each object, we calculate the six SB profiles and their re- 
spective observed luminosities L Q bs(> SB t h), which yields 



a total of 270 data points (LAEs) at each SB t h = 5 x 
10" 20 , 10" 20 , 10" 21 and 10" 25 erg s -1 cm" 2 arcsec -2 . We 
consider the six sides of a simulation cube as part of six sep- 
arate observational fields with comoving side of 5 h~ x Mpc 
each. Structure in the IGM leads to differences in the SB 
profiles for the six sides (as seen in Fig. [4|, which in turn 
leads to a scatter in the observed LFs calculated from the 
different fields. Therefore, LFs are calculated separately for 
each of the observational fields and then averaged in each lu- 
minosity bin to obtain the mean luminosity functions, which 
are shown in Figure [22] 

First we focus on Xi on = 0, to understand the ef- 
fects of SB thresholds on detections (as seen in Figure |9j), 
and determine the luminosity functions. At SB t h = 5 x 
10 -20 erg s -1 cm -2 arcsec -2 only one object is detected 
with an observed luminosity of L Q bs(> SBth) = 3 x 10 38 
erg s -1 , while SB t h = 10 -20 erg s -1 cm -2 arcsec -2 leads 
to the detection of 5 LAEs with observed luminosities larger 
than in the previous case. The LF is thus shifted to higher 
values in both axis. Going one order of magnitude deeper 
in SB leads to a huge increase in the number of detections 
(~ 130 LAEs). This gives rise to a shift in the LF of about 
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two orders of magnitude in both luminosity and number 
density. At SBth = 10 -25 erg s -1 cm -2 arcsec -2 all the flux 
in all the LAEs is observed, returning the intrinsic LF of 
the simulated LAE sample. It should be kept in mind that 
the value of SBth needed to observe the full flux depends 
on the intrinsic luminosity of the sources and neutral IGM 
structure around it. 

In the real universe, both the mean ionization fraction 
and the intrinsic properties of galaxies change. This leads 
to a complex evolution of the LF of LAEs with redshift. 
In addition, the different SB t h from real observations add 
another level of complexity to the redshift evolution of ob- 
served LFs. To help in constraining Xi on with observations of 
LFs, here we try to understand how the LF depends on the 
ionization fraction of the IGM for different choices of SBth- 
For SBth = 10 -25 erg s -1 cm -2 arcsec -2 , the luminosity 
functions for different values of Xi on are equivalent, because 
all the flux from the sources is detected. For lower SB cuts 
instead, the number of detections changes with Xi on , leading 
to an increasing number density of LAEs at a given luminos- 
ity for higher Xi on - The difference between LFs at different 
xion is higher for smaller values of SB t h- It should be noted 
that shallow detections (e.g. 10 -20 erg s -1 cm -2 arcsec -2 ) at 
higher Xi on (e.g. Xi on = 0.89) have similar LFs as deeper but 
incomplete detections (e.g. 10 -21 erg s -1 cm -2 arcsec -2 ) at 
lower Xion (e.g. Xi on = 0). Thus, when comparing luminos- 
ity functions at different redshifts, the SB threshold of the 
observations should be taken into account to obtain a realis- 
tic estimate of the decrease in the number density of LAEs. 
Another point to note is that discrepancies between LFs for 
different values of SBth are smaller at higher ionization frac- 
tions (e.g. x^n = 0.89). Thus only mild differences in LFs 
for orders of magnitude in SB t h at a specific redshift can be 
used as an evidence for high Xi on - At high Xi on , absorption 



outside the edge of the cube (as discussed in Section 6.3.1 ), 
which has not been included in this test, would lead to an 
even larger gap between LFs for different SBth- Neverthe- 
less one has to keep in mind that velocity shifts discussed in 



Section 6.2 (e.g. due to outflows from the galaxy) could also 
lead to milder changes in luminosity functions with changing 
SBth. 

Due to the difficulties in obtaining deep observations 
of LAEs (refer to Section [5]) , we look into the possibility of 
gathering additional information from stacked profiles. This 
was recently done by Steidel et al. (2011) at z = 3 to es- 
timate the low SB extended emission present in individual 



non-detections. In addition, also Zheng et al. (2011) looked 
at estimating total Lya luminosities from stacked profiles 
for different mass and luminosity bins. In our study, we in- 
vestigate the shape of the stacked SB profiles without bin- 
ning for halo properties, to understand how the profiles are 
changed by Xi on - The data set used is the same one discussed 
above. In general, for Xi on = the individual profiles show 
a flatter SB profile with substantial structure due to scat- 
tering through a mostly neutral but structured IGM. For 
higher initial ionization fractions, the strength of the point 
source increases as more photons escape with less scattering 
resulting in a steeper SB profile. These characteristics can 
be exploited to study the mean ionization fraction of the 
universe from stacked LAE SB profiles. 

Figure [23] shows the stacked SB maps from the 270 sim- 
ulated maps of the 45 galaxies described previously. The top 



(bottom) panels refer to the mean (median) values in each 
pixel of the stack. The mean value is generally affected by the 
outliers when the distribution is not gaussian. In our stacks, 
the higher mass haloes act as outliers being brighter and 
more point-source-like due to larger ionized bubbles around 
them, but are fewer in number compared to the majority of 
SB maps of low mass haloes. This leads to the mean profile 
being steeper and brighter than the median profile at each 
Xion- But it is important to note that in both mean and me- 
dian stacks, the trend of the stacked SB profiles with Xi on is 
the same. 

As we have seen in Section [6.3. 1| the higher Xi on is, the 
steeper is the SB profile of a LAE. This trend is more im- 
portant for haloes of lower masses, as the ionization bubble 
around these objects are smaller due to the lower ionizing 
photon flux, leading to a more diffuse SB profile at low Xi on - 
As x^n increases, a larger fraction of the photons from these 
haloes escapes unscattered, leading to a more concentrated 
SB profile. For high mass haloes, even at low Xi on = 0, the 
photons escape with lesser scatter due to larger ionized re- 
gions around them, leading to a steeper SB profile. Due to 
the shape of the halo mass function though, there is a much 
larger number of low mass haloes in a comoving volume of 
the universe. Therefore, when stacking SB profiles of LAEs 
in a volume, the average of the stack is dominated by the 
low mass haloes. Thus the stacked SB profile at low Xi on = 
is flat. The stacked profile has no structural details, because 
the profiles in individual SB maps have a random distribu- 
tion, thus averaging out in the stack. At ionization fractions 
of x^n = 0.5, 0.89, the photons escape easier in all directions 
leading to a steeper SB profile in all directions. This also ap- 
pears when stacking SB profiles, with the stacked SB profile 
becoming; more concentrated with increasing Xion • Thus this 
method of stacking can be used to estimate the steepness of 
the profile and in turn the volume averaged ionization frac- 
tion of the universe at that redshift. As already discussed, 
the velocity shifts could also lead to the steepening of the 
individual SB profile even at Xi on = 0, thus reducing the 
difference in concentration between stacked profiles for dif- 
ferent Xion- But absorption due to IGM outside the cube at 
x^n = would lead to flattening of individual SB profiles. 
This would help in differentiating between Xi on in stacked 
SB profiles. 

In short, a lot of information about the galaxies and 
their surrounding IGM can be obtained from the deep ob- 
servations of LAE SB profiles. A more detailed modelling 
and comparison with observations is the work of a future 
study. 



8 DISCUSSION AND CONCLUSIONS 

The Epoch of Reionization (EoR) is an interesting and im- 
portant event in the history of the universe, the time be- 
ing when H in the gas changed from a mainly neutral to a 
highly ionized state. The details of the process are still un- 
clear and are the topic of a large number of theoretical and 
observational studies. Lyman Alpha Emitters (LAEs), due 
to their strong Lya emission line which scatters for even tiny 
amounts of neutral hydrogen, are used as one of the tools 
to probe the EoR. Modelling LAEs is a complex task with 
a large number of parameters and scales being important. 
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Figure 23. Stacked SB profiles for initial volume averaged ionization fractions of Xi on =0 (left panels), 0.5 (middle) and 0.89 (right). 
The top (bottom) panels refer to the mean (median) values in each pixel of the stacked SB profiles. 



Different studies either focus on a detailed understanding of 
individual aspects of LAEs or use numerical/semi-analytic 
approaches to model the LAEs and compare them to ob- 
servations. In this paper, we focus on the Inter-Galactic 
Medium (IGM) close to the ionizing/Lya source and its ef- 
fect on the observability of LAEs due to Lya RT. The IGM 
close to the source is important: the density, temperature 
and ionization structure controls the Lya radiative transfer 
(RT), determining the intensity and surface brightness (SB) 
distribution along different lines-of-sight. 

To achieve this, we simulated a sample of more than 100 
LAEs, using galaxies and surrounding IGM extracted from 
simulation boxes of 5-30 /i -1 Mpc at z = 7.7. Coupled RT of 
ionizing and Lya photons was performed using CRASHa, 
to determine the ionization structure carved in the IGM by 
the ionizing photons exiting from the galaxies and the spec- 
trum of the Lya photons scattering through the remaining 
neutral hydrogen along different lines-of-sight. The outputs 
were also used to produce Lya SB maps for the lines-of- 
sight perpendicular to the six sides of each simulation cube. 
Analyses were done to study individual objects as well as 
statistical trends. A parameter study was also undertaken 
to understand the dependence of the results on the differ- 
ent, currently uncertain factors involved in the simulations. 

Our main results are: 

(i) Inhomogeneities in the IGM affects Lya RT, lead- 
ing to structure in the simulated SB maps of LAEs. The 
anisotropic escape of Lya photons through the IGM causes 
large variations in the total flux in SB maps along differ- 
ent lines-of-sight for the same object. This leads to ~ 30% 
more scatter in the observed luminosity-mass relation than 



the intrinsic one for our sample of simulated LAEs. This 
also leads to a spread in the values for the escape fraction of 
Lya photons from the IGM: / eS c,Lya,iGM = 0.73 ±0.18 (lcr). 
Note that some lines-of-sight, especially through voids, can 
have /esc, Lya, igm > 1 due to Lya photon contributions from 
other lines-of-sight by scattering and higher probability of 
Lya photon escape towards the observer through voids. 

(ii) Observational campaigns have SB thresholds, SB t h, 
which limit the fraction of the flux observed from the SB 
distribution. At SB t h (e.g. 10 -20 erg s -1 cm -2 arcsec -2 ), 
the observability of the LAE strongly depends on the IGM 
ionization structure and velocity field around the objects. 
Therefore, in observational campaigns at very low SBth, a 
single detection within the observed comoving volume needs 
not be that of the most massive LAE in the region but could 
be a lower mass object seen through a biased line-of-sight. 
At even fainter SB t h, more detections are possible, but the 
observed luminosity values can vary by orders of magnitude 
for a single object depending on the line-of-sight. To ob- 
tain the total flux from the object along each line-of-sight, 
one needs unrealistically deep observations at SBth ~ 10 -25 
erg s -1 cm -2 arcsec -2 . Thus we find that the impact of SB 
thresholds on estimates of observed luminosity is very im- 
portant and needs to be taken into account in the calcu- 
lation of luminosity functions of LAEs from observational 
campaigns (also see Zheng et al. 2010a). 

(iii) One of the main factors which affect the Lya RT 
through the IGM is the wavelength of the Lya photons when 
they escape from the Inter-Stellar Medium (ISM) into the 
IGM. Outflows/winds in the ISM would redshift the pho- 
tons reducing their probability of scattering in the IGM and 
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aiding the escape of Lya photons towards the observer ( Di- 
jkstra & Wyithe||2010L We find that the higher the red- 



shifting of the Lya photons before entering the IGM is, the 
more concentrated is the SB profile, making the detection 
easier at low SB t h and improving the observability of LAEs. 
Therefore, outflows from the galaxies could weaken the im- 
print of reionization on the statistics of LAEs than thought 
previously. But due to the uncertainty in the link between 
the galaxy and outflows, this would add a challenge to esti- 
mating the ionization fraction of the IGM at high redshifts. 

(iv) The ionization structure of the IGM also plays a very 
important role in Lya RT thus determining the SB profiles. 
We estimate the effects of an ionizing background using RT 
simulations with an initial non-zero uniform ionization level. 
At low levels of mean ionization in the IGM (xi on ~ 0), the 
Lya photons undergo scattered diffusion through the IGM 
and the SB profiles thus produced are more extended and 
faint, making it harder to observe the total flux from the ob- 
ject. At high levels of mean IGM ionization (xi on > 0.5), the 
photons scatter less, leading to a more point source-like SB 
profile, making it easier to detect more flux for low SB t h- 
Clustering of sources also affects the ionization structure 
around the object by making the region more ionized, but 
with a complex structure for the ionized bubble due to the 
distribution of neighbouring ionizing sources and IGM gas 
distribution. The qualitative effect on the SB profiles due 
to clustering is similar to that of a uniform ionization case 
with the same mean ionization level. But the detailed dis- 
tribution of flux along different lines-of-sight varies due to 
the differences in the IGM ionization structure. Thus, this 
shows that proper treatment of the ionization structure in 
the region is important for better modelling LAEs. 

(v) Due to the difficulty in achieving very high SBth in 
observational campaigns, stacking of SB maps can be used 
to extract more information from the current LAE samples 
mapped upto relatively lower SBth (e.g. Stei del et al.|201~ 



Zheng et al. 2011). We find that the mean/median stacked 



SB profile of LAEs becomes steeper at higher mean IGM 
ionization levels thus giving an additional technique to use 
LAEs to study EoR. But outflow and halo properties among 
others could also lead to steepening of SB profiles and should 
to be considered when using this technique. 
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